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The dynamical properties of the "Site Frustrated Percolation" model are investigated and com- 
pared with those of glass forming liquids. When the density of the particles on the lattice becomes 
O . high enough, the dynamics of the model becomes very slow, due to geometrical constraints, and 

rearrangement on large scales is needed to allow relaxation. The autocorrelation functions, the 
specific volume for different cooling rates, and the mean square displacement are evaluated, and are 
, found to exhibit glassy behaviour. 

o 
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I. INTRODUCTION 

Frustration plays a central role in many complex systems, such as spin glasses [[l] and glass forming liquids In 
spin glasses frustration arises because ferromagnetic and antiferromagnetic interactions are distributed in such a way 
that the spins cannot simultaneously satisfy all interactions. In glass forming liquids frustration arises when local 
arrangements of molecules kinetically prevent all the molecules from reaching ordered close-packed configurations. 

In this paper we study a lattice-gas model which contains frustration as an essential ingredient. This model is the 
site version of the bond frustrated percolation model, which was defined in the context of the generalization of the 
Kasteleyn-Fortuin Q Coniglio-Klein Q] cluster formulation to the spin glass model Q . 
More precisely consider the Ising spin glass Hamiltonian 
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H = -jJ2<*iS<S j , (1) 



where etj are quenched random interactions which assume the values ±1 with equal probability. By introducing bond 
variables between each nearest- neighbour pair of spins, it is possible to show that the spin glass partition function 
I can be expressed as a sum over bond configurations on the lattice || 

O ' 

Z = Y / ^^ n{C) q N(C \ (2) 

where q — 2 is the multiplicity of the spins, /3/i = ln(e 9/3J — 1), n(C) and N(C) are respectively the number of bonds 
and the number of clusters in the bond configuration C. The asterisk in (||) means that the sum extends over all 
Q 1 the bond configurations that do not contain a "frustrated loop" . This is defined as a closed path of bonds which 
contains an odd number of antiferromagnetic interactions. In this formalism the spin glass problem is mapped on to 
a geometrical problem, which for T = corresponds to pack the highest number of bonds without closing frustrated 
q loops. 

Q , For general values of q the model has been called "q-Bond Frustrated Percolation " || and the partition function 
(||) can be obtained from an hamiltonian model in which, besides the Ising spins interacting in a spin glass way, there 
are Potts spins in each site, interacting ferromagnetically, with multiplicity s = q/2 

The model exhibits two critical points for each value of q: one at higher temperature T p (q), corresponding to the 
percolation transition, in the same universality class of the ferromagnetic Potts model with multiplicity s = q/2, and 
one at lower temperature Tsg(q , ) i in the same universality class of the Ising spin glass transition. This has been 
verified by renormalization methods || , and numerically H . 

It has also been shown [0J|], that for q ^ 2 the percolation transition corresponds to a singularity in the partition 
function, and thus it corresponds to a real thermodynamic transition. 

Each critical point is characterized by a diverging length, associated with the quantities |?]|| : 

Pij=Ptj+Pij ( 3 ) 

and 
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9ij =Pij~ Pij ( 4 ) 

Here pfj (p~) is the probability that 1) sites i and j are connected by at least one path of bonds, and 2) the phase 
r)ij defined as the product over all the signs e mn along the path connecting i and j is +1 (—1) jUl. 

The length £ p associated with the pair connectedness function pij diverges at the percolation critical point, while 
the length £ associated with gfj diverges at Tsg(?) (the bar represents the average over all possible interaction 
configurations {tij})- For q — 2 the quantity g^ in Eq. || coincides with the spin spin pair correlation function. 

However we can give a geometrical interpretation of the second length £. The second transition, like the quantum 
percolation transition |p"l|| , occurs at a bond density higher than the usual percolation transition, due to the interference 
of paths with different phases. At high density the number of allowed configurations is extremely reduced and most 
of the configurations connecting a given pair of sites i and j will have a common path. We call this path a "quasi 
frozen" path, since in a dynamical sequence of configurations exploring the allowed phase space, this path will be 
present most of the time. In conclusion the major contribution to g^ is due to those configurations in which % and 
j are connected by quasifrozen bonds |^,^|,^| and \~g~ij\ roughly coincides with the probability that i and j belong 
tme quasifrozen cluster, and the length £ associated with g 2 

clusters. 

For q = 1 (P) reduces to the bond frustrated percolation model, in which bonds are randomly distributed on the 
lattice, with the only constraint that configurations of bonds which contains at least one frustrated loop are not 
allowed @. 

Thus the model incorporates in the simplest way the concepts of geometrical frustration and could be applied to 
those systems, such as glass forming liquids, where geometrical frustration due to packing problems plays a central 
role. We note that in this model the disorder is quenched while in glass forming liquids this is not the case. However 
at low temperature or high density the relaxation times are so large that the disorder may be considered frozen in. 

To better describe the particle nature of glassy systems, we have introduced in this paper the site version of the 
frustrated percolation model, in which sites can be full or empty. Full sites can be seen as "particles" , and empty sites 
as vacancies. We leave unchanged the underlying structure of the interactions, and the constraint that no frustrated 
loop can be fully occupied. We consider in this paper only the q — 1 case. 

In sect. H we introduce the model, and in sect. [n| we describe two types of Monte Carlo dynamics (thermalization 
and diffusion dynamics), used to simulate the model. We have studied the static properties of the system in sect. 



to the same quasifrozen cluster, and the length £ associated with g~- roughly represents the linear dimension of these 



IV. In particular we have verified that the percolation transition is in the same universality class of the s = 1/2 
ferromagnetic Potts model. The cooling rate dependence of the specific volume is considered in sect. [v|, while in sect. 



VI we analyse the diffusion dynamics of th e model, evaluating the mean square displacement of the particles and the 



diffusion coefficient. Finally, in sect. VH we evaluate the autocorrelation functions of density fluctuations, both in 
thermalization and diffusion dynamics, which show a behaviour typical of glass forming liquids. In the appendix we 
describe in more detail the Monte Carlo procedure. 



II. THE "SITE FRUSTRATED PERCOLATION" MODEL 



In the site frustrated percolation (SFP) model particles are introduced on the vertices of a regular lattice. We assign 
to each edge of the lattice an interaction = ±1, positive or negative according to a random quenched distribution. 
Like in spin glasses a frustrated loop is a closed path of edges which contains an odd number of negative interactions. 
A particle configuration is forbidden only if it contains a fully occupied frustrated loop |l5| (see Fig. |). Thus we 
assign zero weight to all configurations of particles which contain at least one fully occupied frustrated loop, and to 
all other configurations a weight 

W(C) = e^ n{c \ (5) 

where n(C) is the number of particles in the configuration C, \i is the chemical potential of the particles and (3 = 1 /ksT. 
The partition function of the model is given by 

Z = Y J *e^ n[C \ (6) 
c 

where the sum excludes forbidden configurations. 

Note that there is only one independent parameter []p, which can vary from — oo to oo. Frustration prevents the 
system from reaching maximum density p = 1 for (3p — > oo: when (3p varies from — oo to oo, the density p = (n) 
varies from to p m ax < !• 
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The SFP model is expected to have two critical points, in analogy to the bond case ||7|-g|. Fixing the value of p to 
a positive value, and varying the temperature from high to low values, there is a first critical point at a temperature 
T p , corresponding to the percolation transition, and one at lower temperature Tgc, corresponding to a spin glass 
transition. The percolation transition is expected to be in the same universality class of the ferromagnetic s = 1/2 
ferromagnetic Potts model, while in two dimensions Tsg is expected to occur at T = 0. 

At low temperature the model is expected to have the characteristic features of a spin glass system, such as a rough 
free energy landscape, very long relaxation times due to the high free energy barriers, and many ground states at 
T = 0. 



III. MONTE CARLO DYNAMICS 



We have realised two types of Monte Carlo dynamics to simulate the SFP model. The first, which we call "ther- 
malization dynamics" , proceeds through the following steps: 

• pick up a site at random; 

• if the site is filled by a particle, destroy that particle with probability P_, or leave the site filled with probability 

(i-P-); 

• if the site is empty, leave it empty if a new particle placed in that site would close a frustrated loop otherwise, 
if the particle would not close a frustrated loop, create a new particle in that site with probability P+, or leave 
the site empty with probability (1 — P+). 

Looking at the partition function (0), it is easy to verify that the following probabilities of creating or destroying a 
particle satisfy the principle of detailed balance: 

P_ = 1 

p = e /3 M for fin < 0; (7a) 



\p~ = l for fa > 0. (7b) 

This dynamics is clearly ergodic, because we can go from a permitted configuration A to a permitted configuration 
B, first destroying the particles belonging to A, and then creating those belonging to B. 

The difficult step here is to verify that a new particle doesn't close a frustrated loop, since this involves a non local 
check. The procedure is described in the appendix. 

The second type of Monte Carlo dynamics, which we call "diffusion dynamics", starts from a site configuration 
with some density p of particles, obtained by thermalization dynamics, and proceeds by letting the particles diffuse 
conserving their number. At each step a site is chosen at random; if the site is filled, we make an attempt to move 
the particle to a NN site. The move is accepted if the probed site is empty and no frustrated loop is closed. 



IV. PERCOLATION TRANSITION AND EQUILIBRIUM DENSITY 

In this section we analyse the percolation properties of the model. The analysis of the data confirms that the 
percolation transition, for the site frustrated percolation problem, is in the same universality class of the ferromagnetic 
s = 1/2 Potts model, as was formerly verified for the bond problem ||. 

We have used the histogram method for analyzing data |16[ ). For various lattice sizes, we simulated the system 
using thermalization dynamics for 10 temperatures around the percolation point. For each temperature we reached 
equilibrium by taking 10 4 steps, and then evaluated the histograms of the following quantities, taking 10 5 steps: 

• density of particles p; 

• probability of existence of a spanning cluster ; 

• mean cluster size x- 
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The mean cluster size is defined as 



x = ^E s2 ™- ( g ) 

s 

where n s is the number of clusters having size s in the system. 

Using the histogram method |ll|, we evaluated the values of these three quantities for an entire interval of the 
parameter Pp. In Fig. ^ is shown the finite size scaling of and of x, from which it is possible to extract the values 
of the critical temperature, and of the critical exponents v and 7 |17j]. 

The data are perfectly compatible with the values v^ 1 = 0.56, 7 = 3.27, of the s — 1/2 Potts model fisj| . The 
critical (percolation) value for Pp is found to be (Pp) p = 1.2, which corresponds to a density p p ~ 0.6 which is, within 
the error, not different from the percolation density of the standard random site percolation, p p = 0.593 |fl9|| . 

We evaluated the equilibrium density of the SFP, varying the parameter Pp from Pp = —15 to Pp = 15, with 
rate 3 x 10 -5 step -1 . The result is shown in Fig. 0. For comparison we plot also the equilibrium density of the 
unfrustrated model, which is given by p = 1/(1 + e - ^). This curve fits quite well the low density part, since for low 
density the effect of frustration is negligible. The high density part instead is better fitted by the function 



_ Pmax /„n 

9 ~ 1 + e -«-b(/3M) 1 ' 

where p max = 0.788, a — 0.323, b — 0.719. Note from Fig. || that the crossover occurs at a value of Pp, close to the 
percolation threshold. 



V. COOLING RATE DEPENDENCE OF THE DENSITY 



The effect of frustration prevents the system from easily reaching the equilibrium density, especially at very high 
values of the parameter Pp, that is for Pp — > 00. 

Fixing p to a positive value, and cooling the temperature from high to low values with a finite cooling rate, there 
is a temperature Tq at which the system goes out of equilibrium. Lowering further the temperature, the density 
remains constant to the value corresponding to the temperature Tq. 

Fig. ^ shows the temperature dependence of the density for various cooling rates. Temperature goes from ksT/p — 
0.5 to ksT/fi = 0, and cooling rates range from ksT/p = 10 -2 step -1 to 10 -7 step -1 . One step is one update per 
site. 

This behaviour is experimentally observed in supercooled glass-forming liquids, at the calorimetric glass transition 
temperature Tq, when structural relaxation times become greater than the experimental observation times pp| . As 
in glass forming liquids, in the model we observe that the faster the cooling rate, the greater are the glass transition 
temperature and the specific volume. 



VI. THE DIFFUSION DYNAMICS 



We have then studied the diffusion dynamics (conserved number of particles) of the system. The diffusion process 
is severely hindered by kinetic constraints at high density, and particles can diffuse through the system only by a 
large scale cooperative rearrangement of many particles. At low density we expect that the effect of frustration on 
the diffusion is very weak because of abundance of holes, and particles diffuse freely on the lattice. 

We have evaluated the mean square displacement (MSD) (Ar(t) 2 ) as a function of time t, as shown in Fig. |[ Each 
curve is obtained averaging over all the particles and over a time interval of 5000 steps. As expected for low density 
the curves show a linear behaviour, which corresponds to normal diffusion. For high densities, the MSD reaches a 
plateau and then becomes again linearly dependent on time. This behaviour is typical of glass forming liquids and is 
also observed in molecular dynamics simulations |^],^2j . The crossover from the normal behaviour to the anomalous 
diffusion in our model again occurs at a density value close to percolation. 

We have extrapolated the diffusion coefficient values D — (Ar(t) 2 )/t from the long time regime which are shown 
in Fig. H The extrapolated values are very well fitted by a power law as function of density. 

d ~ ( P - Pc y (10) 

with p c ~ 0.7874 and 7 ~ 1.5. 
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Note that p c is equal within the errors to p max , the density corresponding to the ground state, as extracted from 
the fit (^) . Thus the dynamic singularity coincides in this model to the spin glass transition, which in two dimensions 
occurs at T = 0. If we set p c = p max , using Eq. (S) we obtain from Eq. ( |l(i| ) 

D~(l - ^1 ~ e - b T(^) (11) 

which shows that the diffusion coefficient goes to zero with an Arrhenius law for T — -> 0. 

The picture which emerges from the model is that particles are trapped most of the time in cages formed by their 
NN, and can diffuse only into localized pockets. Only after a very long interval of time they have the opportunity of 
escaping from the cage, falling in a "next neighbour" cage. This process forms a kind of random walk, reflected in 
the long time linear dependence of the MSD functions. 

This picture is enforced by direct observation of the particles moving on the lattice. We have taken "snapshots" of 
the system at particular times, for diffusion dynamics at two high values of the density. Some particles leave tracks 
on the lattice, so we can analyze also the path they have walked. In Fig. [?] is shown the snapshot at p — 0.783 after 
2000 steps (a), after 4000 steps (b)and at p = 0.785 after 2000 steps (c) 

We mark with dots the particles that moved at least one time (non frozen domains) during the time of the simulation; 
the filled circles represents frozen domains. Clusters of frozen particles prevent other particles from moving freely. 
We note that frozen particles are clustered and we can distinguish liquid-like zones (holes and mobile particles) from 
solid- like zones (frozen clusters). Frozen clusters are eroded as the time of the simulation gets longer, while liquid-like 
zones grow. 

At the very high density, frozen particles remain blocked for extremely long times. The frozen cluster in Fig. [7](c), 
corresponding to density p — 0.785, is still frozen after 2 x I0 5 steps which corresponds to the maximum time we have 
observed. We have checked that the large frozen cluster does not percolate in either direction, even if the density 
of the frozen particles is much greater than the percolation density of the random site percolation. We expect that 
the frozen particles in fact percolate at T — ( p = p ma x) since the linear dimension of the frozen cluster should 
correspond roughly to the correlation length associated to (^), which for d = 2 is expected to diverge at Tg G = 0, 
as explained in the introduction. 



VII. THE RELAXATION FUNCTIONS 

An important property that characterizes glassy behaviour is the form of the relaxation functions pq ]. We have 
evaluated the relaxation functions of the system, in the thermalization and in the diffusion dynamics. For each 
temperature, we reach the equilibrium, and then evaluate the relaxation functions averaging on a time interval of 
1 3 - I0 6 steps. 

Fig. H shows the density-density autocorrelation function in the thermalization dynamics, defined as 

_ (Sp(t)Sp(0)) 

m - W) ' ( } 

where Sp(t) — p{t) — (p). 

For diffusion dynamics, we have studied the autocorrelation function of the density fluctuations 

^fe(l) - 77 J- v , (13) 

where 

a 

and r a is the position of the a-th particle in units of lattice constant. The wave vector can take the values k — ^-n, 
where n has integer components n x and n y ranging in 0...L/2. The autocorrelation functions corresponding to 
n x = L /2 and n y = are reported in Fig. 



Fig. |10| shows the self-part of the relaxation function in the diffusion dynamics, defined as 

^ s w- 1 (E e * (rQ(t) " r ° (0)) )' ( 15 ) 
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where n is the number of particles. 

At low temperatures, below the percolation threshold, we observe the onset of non-exponential decay. This is usually 
the sign of non-stochastic cooperative relaxation in the system. We emphasize that, for our model, the percolation 
transition corresponds to a real thermodynamic transition, and therefore it is possible to expect, below this point, a 
change in the dynamical properties of the system. 

At very high density it is evident the existence of different time regimes as predicted by the mode mode coupling 
theory p^ , p^] of supercooled liquids and observed both in some molecular dynamic simulations and in some experi- 
mental measurement |23f| on glass forming liquids. There is a first short time relaxation, corresponding to relaxation 
inside non frozen domains surrounded by a frozen "cage" , and a long time regime (a relaxation) , corresponding to 
structural rearrangement and final decay to equilibrium. 



VIII. CONCLUSIONS 



We have studied a frustrated lattice gas model which it has been called site frustrated percolation since percolation 
plays an important role. This model, despite its simplicity, shows a complex dynamical behaviour. Disorder and 
frustration are its basic ingredients. The constraint of frustration prevents the system from easily reaching high 
density, and inhibits the motion of particles at high density. 

We observed that the dependence of the volume as function of the temperature varies strongly with the cooling rate, 
qualitatively in the same way as observed in real glass forming liquids. The motion of particles in diffusion dynamics 
is severely hindered by the frustration constraint, and the relaxation process can take place only through large scale 
rearrangement. The cooperative nature of the motion is reflected in the behaviour of the mean square displacement 
and in the relaxation functions of the system. The plateau observed in the self part of relaxation functions (Fig. |Io| ) 
and in the mean square displacement (Fig. ||), are intrinsically connected, and are the sign of the glassy nature of the 
dynamical properties of the model. In all these phenomena the crossover from a normal behaviour to an anomalous 
behaviour occurs close to or at the percolation threshold. 

Since SFP seems to well describe the glass transition in glass forming liquids, the model suggests that the presence 
of a percolation type transition may be a general feature below which frustration effects start to be manifested. This 
transition may be responsible of various precursor phenomena |l3| |, such as onset of stretched exponentials, breakdown 
of Stokes-Einstein relation, presence of spatial heterogeneity. 

The presence of a percolation transition well above the glass transition has recently been discovered by Tomida and 
Egami [p8| in a molecular dynamic simulation of monatomic liquids. It is also interesting to note that Kivelson et al. 
p9[ showed that the viscosity of 15 glass forming liquids could be collapsed on one single curve, by assuming only one 
characteristic temperature well above the glass transition. 

Finally the ideal glass transition temperature, characterized by the divergence of the inverse of the diffusion coeffi- 
cient, numerically seems consistent with the divergence of the static length £ associated to <?y ([|). 

Although we have discussed the site frustrated percolation model in the context of glassy systems, the model is 
rather generic and may be applied to other systems, where geometrical frustration plays an essential role. In fact 
recently the model has been successfully applied also to granular materials fl30f . 

We would like to thank V. Cataudella, F. di Liberto, S.C.Glotzer, M. Nicodemi and U. Pezzella for interesting 
discussions. This work has been supported in part by CNR. 



APPENDIX 



As mentioned in sect. Ill, we describe here the procedure to check whether a new particle, added to the system in 
a given configuration, closes or not a frustrated loop. 

Every nearest neighbour (NN) particle (occupied site) of the empty site, that has been probed to be filled, belongs 
to a (not necessarily distinct) cluster of connected sites. If two particles NN to the empty site belong to the same 
connected cluster, then a new particle filling that site closes a loop. More precisely, if z n is the number of particles 
NN to the empty site, and z c is the number of distinct clusters to which they belong, than the total number of new 
loops closed by the new particle is A = z n — z c . 

The algorithm proceeds as follows: 

• we count the number z n of particles NN to the empty site which we want to fill, and mark each of them as the 
root of a distinct cluster, so at the beginning z c = z n , and A = 0; 
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• we grow every cluster in parallel, adding to it in recursive way its still not visited NN particles, and marking 
them as belonging to that cluster; 

• if no NN particle to a cluster is found, that cluster is marked as "burnt" , and is not considered anymore; 

• if two clusters collide, we say we have found a loop, and the two clusters merge to form a single cluster; then A 
is increased by one, and z c is decreased by one. 

Every new visited site is marked in two ways: with the label of the cluster to which it belongs, and with the "parity", 
that is the number of antiferromagnetic interactions walked through starting from the initial empty site. We stop the 
iteration when one of the following circumstances happens: 

• two clusters collide and the parities don't correspond (that is one is odd and the other is even): in this case we 
have found a frustrated loop; 

• the number of burnt clusters equals z c — 1 (only one cluster is non-burnt): in this case no other loop can be 
found. 

When the density of particles is low, this algorithm is very fast, since the clusters are very low sized, and are burnt 
within few iterations. On the other side, when the density is high, the maximum number of independent loop is 
found quickly as well. The algorithm can become slow when the density of particles is intermediate, notably near the 
percolation transition, because clusters in this case are very ramified. In Fig. |ll| the CPU time needed to do a single 
site update is shown for a square bidimensional system, and for various lattice dimensions. The maximum time is 
reached near the percolation transition, and scales as 

woe TV -™ 5 , (16) 

where TV = L 2 is the total number of sites. 
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FIG. 1. (a) Site configuration that does not contain frustrated loops (permitted); (b) site configuration that contains a 
frustrated loop (forbidden). 



FIG. 2. (a) Finite size scaling of the probability that a spanning cluster exists; (b) finite size scaling of the mean cluster size. 

FIG. 3. The equilibrium density of the SFP as function of f3fi (circles), together with the fit function (lower curve). Upper 
curve is the density of the unfrustrated model. 



FIG. 4. Temperature dependence of the specific volume in the SFP model (32 x 32 square lattice), for various cooling rates. 
From upper curves to lower ones: k B f/fi = 10~ 2 , 10~ 3 , 10" 4 , 10~ 5 , 10~ 6 , 10~ 7 step" 1 . 



FIG. 5. Mean square displacement in the SFP model on a square lattice with L = 32, for densities (from upper curves to 
lower ones): p = 0.452, 0.586, 0.730, 0.756, 0.766, 0.777, 0.784, 0.785. 



FIG. 6. Long time diffusion coefficients D as function of density. The fit function is D = 2.14 x \p — 0.7874| 153 . 

FIG. 7. Path of tagged particles at density p = 0.783 after 2000 steps (a), after 4000 steps (b), and at p = 0.785 after 2000 
steps (c). Dots are particles that have moved at least one time in the time indicated (non frozen domains); the filled circles 
represent particles which have never moved in the time indicated (frozen domains) The frozen domains in (c), are still frozen 
after 2 x 10 5 steps. 

FIG. 8. Relaxation functions of the site density, for the thermalization dynamics, on a L = 36 square lattice, for temperatures 
between e 13 ^ = 1 and e 13 ^ = 1000, corresponding to equilibrium density (from lower to higher relaxation times): p — 0.450, 
0.639, 0.715, 0.782, 0.784. 

FIG. 9. Relaxation functions of the density fluctuations, for the diffusion dynamics, for densities p = 0.452, 0.639, 0.777, 
0.784, on a L = 36 square lattice and k x — tt, k y — 0. 



FIG. 10. Self part of relaxation functions of the density fluctuations, for the diffusion dynamics, for densities p = 0.452, 
0.639, 0.757, 0.783 on a L = 36 square lattice, with k x = tt, k y = 0. 



FIG. 11. CPU times needed to make a single site update for lattice sizes L = 24, 32, 40, 48, 56, 64. 
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